function [loglik] = mykalman(data,T,Z,Q,H,xini,P,alf)
%Compute kalman filter
% state(t) = T*state(t-1) + err1(t), Eerr1(t)*err1(t)'=Q
% data(t) = Z*state(t) + err2(t), Eerr2(t)*err2(t)'=H

%Written by Gianni Lombardo. Modified by Roberto Motto

bigt = size(data,2);
n    = size(data,1);

% initialize likelihood with the Gaussian constant
loglik = (n*bigt*0.5)*log(2*pi);% negative of costant part of likelihood times #series*lenght

scfc=1;

for i=1:bigt
    
    Ft = Z*P*Z'+ H; % observables' forecast error
    
    if rcond(Ft) > 1e-8;
        invFt = inv(Ft);
    else
        loglik = -10e12;
        return
    end
    
    ut = data(:,i)-(Z*xini) - alf;
    
    loglik = loglik + (0.5)*(log(abs(det(Ft))) + ut'*invFt*ut);  %likelihood function
    
    % start the filter
    Kt   = T*P*Z'*invFt; %Kalman gain (see Hamilton p. 381)
    xini = T*xini + Kt*ut; 
    %L    = T - Kt*Z;
    %P    = T*P*L'+ Q; %update of MSE of projection of state variables
    P = T*P*T' - Kt*Z*P*T' + Q;
end
loglik = -loglik;